%% simulations over lifcycle
i_exit=zeros(n_sim,T);
n_exit=zeros(n_sim,1);
b_exit=zeros(n_sim,1);
z_exit=zeros(n_sim,1);
t_exit=ones(n_sim,1)*T*dt;
inde_z_hist=zeros(n_sim,T);
inde_b_hist=zeros(n_sim,T);
z_hist=zeros(n_sim,T);
b_hist=zeros(n_sim,T);
l_b_hist=zeros(n_sim,T);
f_hist=zeros(n_sim,T);

obse=0;fund_count=0;sim_count=0;event_study = 1000*ones(n_sim,T);
for ii=1:n_sim
    bd=rand;
    inc_draw=rand;
    if inc_draw>incid
        b0=b_up_unfund*(bd<=q_hig)+b_md_unfund*(bd>q_hig)*(bd<=q_hig+q_mid)+b_dn_unfund*(bd>q_hig+q_mid);
    else
        b0=b_up_fund*(bd<=q_hig)+b_md_fund*(bd>q_hig)*(bd<=q_hig+q_mid)+b_dn_fund*(bd>q_hig+q_mid);
    end

    z0=zlow*(bd<=q_hig)+zmid*(bd>q_hig)*(bd<=q_hig+q_mid)+zhig*(bd>q_hig+q_mid);

    t=0;inde_alive=1;t_sim=0;
    while t_sim<T && inde_alive==1
        t=t+dt;t_sim=t_sim+1;

        ell_t       = min(100,fnval(cub_spl_l,[b0;A]));
        drift_I_t   = min(0.12,max(0,fnval(cub_spl_d,[b0;A])));
        drift_dlb   = dt*(ell_t/b0+sigma^2/2-rho-drift_I_t);
        drift_dlz   = dt*(drift_I_t-sigma^2/2);
        dw  = randn;
        l_b = log(b0) + drift_dlb-sqrt(dt)*sigma*dw;b0  = exp(l_b);
        l_z = log(z0) + drift_dlz+sqrt(dt)*sigma*dw;z0  = exp(l_z);
        if b0<b_bar
            b_hist(ii,t_sim)  =b0;
            l_b_hist(ii,t_sim)=log(b0);
            z_hist(ii,t_sim)  =z0;
        else
            neg = rand;
            if neg >phi || b0*a>b_bar
                inde_alive= 0;
                t_exit(ii)= t;
                i_exit(ii,t_sim)= 1;
            else
                b0=b0*a;
                b_hist(ii,t_sim)  =b0;
                l_b_hist(ii,t_sim)=log(b0);
                z_hist(ii,t_sim)  =z0;
            end
        end
    end
end

%% yearly averages
tot_years       = T*dt;
mea_exit_surv_y = zeros(tot_years,1);
mea_exit_ent_y  = zeros(tot_years,1);
b_surv          = zeros(tot_years,1);
z_surv          = zeros(tot_years,1);
b_surv_ent      = zeros(tot_years,1);
surv_t = zeros(tot_years,1);

n_surv=n_sim;inde_surviv=ones(n_sim,1);
for tt=1:tot_years
    t0= (tt-1)/dt+1;
    t1= tt/dt;
    n_died        = sum(sum(i_exit(inde_surviv>0,t0:t1)));
    exit_surv(tt) = n_died/n_surv;
    b_surv(tt)    = sum(sum(b_hist(inde_surviv>0,t0:t1).*z_hist(inde_surviv>0,t0:t1)))/sum(sum(z_hist(inde_surviv>0,t0:t1)));
    z_surv(tt)    = mean(mean(z_hist(inde_surviv>0,t0:t1)));
    n_surv        = sum(1-i_exit(inde_surviv>0,t0));
    surv_t(tt)    = n_surv/n_sim;
    inde_surviv   = b_hist(:,t1);
end
z_all     = mean(z_hist(:,end));

exit_prob_yearly  = exit_surv+delta*(1-exit_surv);
